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1.  Introduction 

A  triangular  matrix  reveals  its  eigenvalues  along  the  main  diagonal. 
By  Schur's  lemma  any  square  complex  matrix  is  unitarily  similar  to 
an  upper  triangular  matrix  with  the  eigenvalues  arranged  in  any 
desired  order  along  the  main  diagonal.  It  follows  that  any  real  square 
matrix  is  othogonally  similar  to  a  real  block  upper  triangular  matrix 
in  which  each  2x2  block  on  the  diagonal  corresponds  to  a  pair  of 
complex  conjugate  eigenvalues.  The  Householder- QR  algorithm  is  a 
stable,  efficient  algorithm  that  produces  a  Schur  form.  However  the 
ordering  of  the  eigenvalues  that  the  QR  algorithm  produces  may  not 
be  suitable  for  certain  purposes,  such  as  computing  the  exponential  of 
the  original  matrix.  There  are  programs  in  the  library  EISPACK  that 
compute  this  real  Schur  form.  See  section  23J5  of  |E3S,1976]. 

This  investigation  presents  and  compares  all  the  attractive  methods 
we  can  think  of  for  performing  orthogonal  similarity  transformations 
that  preserve  block  triangular  form  but  rearrange  the  eigenvalues. 
This  is  a  fairly  straightforward  task  but  it  is  always  a  challenge  to 
try  and  keep  down -three  conflicting  costs:  round  off  error,  execution 
time,  and  program  length. 

We  give  some  attention  to  the  task  of  swapping  adjacent  diagonal 
blocks  of  orders  p  and  q  but  our  main  concern  is  with  the  case  p=q=2. 
We  use  capital  letters  to  denote  matrices.  Fortran  programs  are  given 
at  the  end. 

Before  plunging  into  details  we  describe  the  methods  in  brief  general 
terms.  Algorithm  0  (GW.  Stewart):  Swap  adjacent  blocks  using  one 
or  two  QR  steps  with  a  pre-determined  shift  to  force  the  ordering  of 
the  eigenvalues  of  the  new  blocks.  Algorithm  1:  Swap  adjacent 
blocks  as  needed  using  an  explicit  orthogonal  similarity 
transformation.  At  most  4  rows  and  columns  will  be  modified  at 
each  swap.  Algorithm  2:  Swap  adjacent  blocks  using  Householder 
transformations.  For  swapping  a  pair  of  2x2  blocks  two  Householder 
transformations  are  needed. 

The  table  below  compares  the  algorithms  for  code  length  and  running 
time  THe  remainder  of  that  paper  is  concerned  with  accuracy 
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Algorithm  Fortran  line  Speed  ratio  (The  ratio  was 
number  count  determined  by  runs  on  9x9  matrices) 


0 

165 

12? 

1 

366 

1D0 

2 

301 

1.15 

Table  1 


We  were  encouraged  to  present  our  programs  and  results  by  Dr. 
Sven  Hammarling  (of  NAG,  Inc.,  Oxford)  who  showed  us  work  on 
swapping  that  he  had  begun  in  cooperation  with  Dr.  J.  Dongarra 
(Argonne  National  Laboratory)  and  the  late  Prof.  J.  H.  Wilkinson. 


2.  The  Software  Economiser  (EXCHN6) 

Consider  a  submatrix  of  the  form 

A  B 

1 

0  A 

2 

where  A*  and  A2  are  2x2  diagonal  blocks. 

Algorithm  Q  (called  EXCHNG  in  [Ste,  19761). 

1.  An  implicit  double  shift  is  determined  from  A*. 

2.  An  arbitrary  QR  step  is  performed  to  destroy  the  triangular 
form  and  put  the  matrix  into  Hessen  berg  form. 

3.  A  sequence  of  double  QR  steps  using  the  shift  from  step  1.  The 
The  eigenvalues  of  the  first  block  will  emerge  in  the  lower 
part  of  the  array  occupied  by  both  blocks,  usually  in  one  or 
two  steps 

Remark  1.  The  alogrithm  discards  the  information  that  there  are 
two  pairs  of  conjugate  complex  eigenvalues.  Stewart  modifies  the 
standard  QR  program  so  that  a  supplied  initial  shift  may  overwrite 
the  usual  Francis  shift  at  the  first  step.  Such  an  algorithm  would 
converge  in  one  or  two  steps. 
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3.  General  Theory 

Consider  the  block  upper  triangular  matrix 

A  is  pxp, 


A  B 

1 

0  A 


'  A2  is  q*q. 


(1) 


Throughout  this  paper  we  assume  that  and  A2  have  no  eigenvalue 
in  common.  It  follows  that  there  exists  a  unique  pxq  matrix  X  such 
that 

AjX  -  X  A2  =  B.  (2) 

This  is  called  Sylvester's  equation.  It  follows  that 

A 


A  B 

1 


Ip  -X 

0  I. 


-*  *p 

«q° 


(3) 


DEFINITION.  An  orthogonal  (p+q)x(p+q)  matrix  H  is  said  to  swap  A* 
and  A2  if 


H 


A  B 

1 


H 


A  B 

2 


(4) 


where  A*  is  similar  to  Aj,  i=l,2. 


Lemma  1.  An  orthogonal  (p+q)x(p+q)  matrix  H  swaps  A^  and  A2  if, 
and  only  if. 


(5) 


-X 

M 

H  • 

V 

- 

2 

0 

for  some  invertible  qxq  M2  where  X  is  defined  in  (2). 
Note  that,  since  H  is  invertible, 

M 

=  q 

Consequently  M2  is  qxq  and  must  be  invertible. 


M 

-X 

rank 

2 

=  rank 

V 

0 

Proof.  If  H  satisfies  (5)  then  for  some  qxp  W  and  pxp  Mj, 

[  -X  I  1  [  M  W  1 


*q  0 


M  W 

2 

0  M 


and,  since  both  matrices  on  the  left  are  invertible  so  are  Mj  and  M2. 
Thus 

A  B  _  -X  In  1  A  0  1  0  I  j  _ 

H  •  1  •  HT  =  H-  p  .  2  .  .hT 

0  A2  Iq  0  0  A1  IpX  °  ' 


M  W 

2 

0  M 


A  0 

2 

0  A 


-1  -1  -1 

M  -M  -W-M 

2  2  1 

0  M^ 


A  B 

=  2 

—  # 

0  A 

1 

where 

Aj  =  Mj-Ai-Mf1 ,  i*  1,2;  B  =  (WAi  -  • 

Conversely,  if  H  swaps  A^  and  A2  then  there  exist  M2,  and  W 
such  that 

I** 6LlM2  wl.h °l.l< 


*2  ® 

m2 

V 

** 

0  A 

lj 

0 

M 

lj 

=  H 


0  A 


-XIAO  01 


2  -1  1 

M 

1 


LOO  A.  1  LX 


q  ht 

*n 


It  follows  that 


M 

W 

-1 

-X  L 

2 

0 

M 

•  H- 

p 

I  0 

lj 

q 

commutes  with 


A  0 

2 

0  A 


Since  Aj  and  A2  have  no  eigenvalues  in  common  D  must  be  a 
polynomial  in  diag(A2,Ai).  See  [Gant.  vol.  1,  p222J. 


-X  I 

M  W 

p 

2 

I  0 

0  M 

q 

lj 

must  be  block  upper  triangular.  This  establishes  the  converse 

QED 


5 


2.  An  othogonal  H  that  swaps  Aj  and  A2  must  have  the 


form 


where 


H  = 


-1 

T 

c  0 

-X  Ia 

2  , 

• 

q 

0  c  1 

1-  X 

1 

T) 

T 

T 

:  -c  = 

1 

♦  X  -X 

2  2 

q 

T 

T 

:  -c  = 

I 

♦  XX 

1  1 

p 

(6) 


(7) 


Proof.  Write  Cj  for  multiply  (5)  by  HT  and  use  the  orthogonality 
of  H  to  find 


-X 

-X 

_T 

_  1 

V 

y 

T 

=  H  -H- 

II 

C2 

n 

• 

0 

1 

0 

Transposing  reveals  the  first  row  of  H.  The  second  row  follows  by 
orthogonality. 

QED 


Remark  1.  There  are  infinitely  many  choices  for  C\  and  0^  that 
satisfy  (7).  See  Section  4.3  for  more  details. 


Remark  2.  One  of  our  implementations  uses  the  form  in  (6)  explicitly. 
The  block  rows  are  orthogonal  by  their  form,  so  it  is  the  accuracy 
with  which  (7)  is  fulfilled  that  determines  the  orthogonality  of  the 
computed  H.  An  alternative  implementation  starts  from  (5)  and 
seeks  H  as  a  product  of  elementary  reflectors  (also  known  as 
Householder  matrices). 


The  key  blocks  of  the  transformed  matrix  can  be  found  explicitly. 
Using  (3),  (6),  and  (7)  it  is  not  difficult  to  see  that 


~  T  -T 

A  =  C  -A  -C 

2  2  2  2 


A  =  C  -A  -C 

1  111 


T  *T  *T  _  4  T 

B  =  C  -A  -X  -C  1  -  C  -X  -A  -C 


2  2 


1  1 


(8) 
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4.  Implementation  details 


4.1.  Standardized  Real  Schur  Form 

The  Schur  form  of  a  matrix  is  not  unique  and  the  real  Schur  form  of 
a  real  matrix  offers  even  more  freedom.  We  urge  the  adoption  of  the 
following  conventions. 

i)  2x2  diagonal  blocks  are  used  exclusively  for  complex 
conjugate  pairs  of  eigenvalues,  not  for  distinct  real 
eigenvalues. 

ii)  The  diagonal  elements  of  2x2  diagonal  blocks  are  mode  equal. 
This  value  is  the  real  part  of  each  eigenvalue. 


Consequently  we  advocate  the  form 

«  t> 

¥  oc 


P*¥  <  0  . 


The  off  diagonal  elements  of  the  2x2  diagonal  blocks  cannot  always  be 
made  equal  in  absolute  value  but  they  must  be  opposite  in  sign.  To 
guar  ant  uniqueness  one  may  require  fr  and  ¥  to  satisfy  0  <  ¥  1 
but  that  is  nol  essential.  Note  that  the  eigenvalues  are  a +  ./&•¥ 

The  use  of  a  standard  real  Schur  form  facilitates  the  swapping  of 
diagonal  blocks  as  well  as  ensuring  that  the  real  parts  of  all 
eigenvalues  are  held  on  the  diagonal  of  the  real  Schur  form. 

If  a  given  real  Schur  form  does  not  have  its  eigenvalues  ordered 
appropriately  down  the  diagonal  then  some  swapping  of  diagonal 
blocks  will  be  needed.  However  the  task  is  considerably  simplified  by 
the  fact  that  no  block  has  order  exceeding  2.  Any  configuration  of 
eigenvalues  can  be  reached  by  swapping  adjacent  diagonal  blocks  and 
this  is  the  task  we  consider  below. 


Here  is  a  method  (cf  section  4.5)  to  put  2x2  diagonal  blocks  into 
standard  form.  Let 


A 


Define  a  reflection  P 
PT  =  P  =  I 


by 

-cos($) 

sin($) 


all  al2 
a21  a22 


sin($) 
cos(a)  J 


>  *  =  2  tan"  ( 


all  -  a22 
al2  ♦  a21 
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Write  c=cos(d)  and  s=sin($).  It  is  not  difficult  to  see  that  PAP 
transforms  A  to  a  standard  form: 

all*a22 
2 

a!2  -  £.411^22 


42  Solving  A,X  -  XA2  =  B 

Considerable  attention  has  been  paid  to  the  general  case  of  this 
equation,  now  known  as  Sylvester  s  equation.  See  [B&S,1972J  and 
{G,N,&vL,  19791 .  When  A^  and  A2  are  either  l*i  or  standardized  2x2 
matrices  the  solution  can  be  given  explicitly  using  stable  formulae 

In  an  earlier  unpublished  report  (Pa, 1977)  we  advocated  scaling  X, 
i.e.,  we  solved 

AjX  -  XA2  =  vB 

with  *  chosen  so  that  HXllssl.  Further  analysis  shows  that  this 
caution  is  unnecessary.  There  is  no  danger  in  working  with  X  of  large 

norm  provided  that  HXll2  does  not  overflow.  Moreover  if  HXll2  does 
overflow  then  the  blocks  should  not  be  swapped  because  a  tiny 
pertubation  will  give  the  new  and  A2  at  least  one  common 
eigenvalue. 

Our  algorithm  for  solving  the  Sylvester  equation  is  called  TXMXT  (for 
TX-XT)  and  is  described  in  Appendix  A,  see  also  the  program 
Appendix  C. 

4.3.  The  Choleski  Factor 

If  the  explicit  orthogonal  H  described  in  Section  3  is  to  be  used  then  it 
is  necessary  to  solve  the  equations  (7)  for  and  C2.  We  can  see  no 
reason  to  avoid  the  Choleski  factorization.  The  formulae  are  given 
below.  When  presenting  the  algorithm  in  detail  we  write  xij  for 
X(i,j).  Recall  equation  (7): 

T  T 

C  -C  =  I  + X  -X  . 

2  2  q 

T  T 

C  -C  =  I  *  x-x  . 

1  1  p 
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Algorithm  1. 

An  H  is  found  explicitly  in  the  form  of  (6)  to  swap  and  A2. 
The  choice  for  and  C2  in  (7)  are  the  Choleski  factors. 


Case  1.  X  is  lxl, 
Case  2.  X  is  1x2, 


then  C2(i,l)  =  <^(1.1)  =  Jl+xll2  . 
then  Ci(l,l)  =  N/l+xll2*xl22,  and 

'  0  •  j 


xll*xl2 


i.(^>2 


l*xll 


Case  3.  X  is  2x1,  then  C2(l,i)  =  v/l*xil2+x212,  and 

*  0 

xll«x21  1 21)2 


l*xll‘ 


Case  4.  X  is  2x2,  let  8  =  xll-x22  -  xl2-x21,  V  -  ^/l+xl  l2*xl22,  and 
k  =  N/l*xll2+x212  ;  we  have 


C  = 
1 


xll*x21+xl2*x22 


and 


V 


xll-xl2+x21*x22 


2  2  2 

x21  *x22  *6 


1  ♦ 


2  2  2 

x!2  *x22  ♦& 


1  ♦ 


4.4.  Representing  H  as  a  product  of  two  reflectors 

The  explicit  form  of  H  in  Section  3  is  not  mandatory.  W.  Kahan 
suggested  using  two  reflections  instead  Here  are  the  details  First 
some  notation:  A  nxn  reflection  (or  Householder  matrix)  can  be 

represented  as  I  -  uu*7d,  where  I  is  the  nxn  identity  matrix,  u  is  a 
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n- vector,  and  d  =  HlluJI2.  We  use  the  fact  that  if  u=x+y  and  11x11= llyll, 
then  (I-uuT/(H||ul|2))’X=-y. 


Algorithm  2. 

An  H  is  found  implicitly  in  the  form  of  either  a  reflector  or  a 
product  of  two  reflectors  to  swap  A*  and  A2.  The  reflectors)  are 
determined  as  follows: 


Case  1.  X  is  lxl.  Let  sx  =  signUlD-^/l+xii2  .  We  seek  a  reflection 
H  so  that 


-xll 

-sx 

1  . 

0  . 

The  special  form  of  H  leads  to 

If  xll/sx  i  05,  then  ul  =  sx  -  xll;  else  ul  =  l/(sx+xll) 
u2  =  1 


d  =  ul-sx 


Case  2.  X  is  1x2.  Let  sx  =  -$ign(xi2)\/l+xll2+xl22  .  We  observe 
that  if  a  reflector  H  satisfies 


1 

0 

H- 

xll 

= 

0 

xl2 

-sx 

then  it  satisfies  (5).  The  proof  is  left  to  the  reader.  The  special  form 
of  H  leads  to 
ul  =  1 
u2  =  xll 

if  -xl2/sx  i  0.5,  then  u3  =  xl2+sx;  else  u3  =  (l+xll2)/(sx-xlZ) 
d  =  u3-sx 

Case  3.  X  is  2x1.  Let  sx  =  signfxlD^l+xll^xZl2  .  From  (5)  we 
seek  a  reflection  H  so  that 


-xll 

-sx 

H- 

-x21 

= 

0 

1 

0  . 
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The  special  form  of  H  leads  to 

If  xll/sx  i  05,  then  ul  =  sx  -  xli;  else  ul  =  (i*x212)/(sx+xli) 
u2  =  -x21 
u3  =  1 
d  =  ul-sx 

Case  4.  X  is  2x2.  Two  reflections  Hi  and  H2  are  required.  In  (5)  let 
M2  be  upper  triangular  and  H  =  H2’Hi.  First  define 

sx  =  sign(x)*v/l+xll2+x212  . 

We  seek  a  reflection  Hi  =  1  -  uuT/d  so  that 


H 


From  the  special  form  of  H  we  have 

If  xll/sx  <  0.5,  then  ul  =  sx  -  xll;  else  ul  =  (l+xZl2)/(sx+xll) 
u2  =  -x21 
u3  =  1 
u4  =  0 
d  =  ul-sx. 

Next,  define  an  intermediate  vector  y  by 


-xll 

-sx 

-x21 

0 

1 

0 

0 

0  . 

yi 

-xl2 

y2 

=  H  • 

-x22 

y3 

1 

0 

y^ 

1  . 

One  can  verify  that 

yl  =  -(xli*xl2+x21-x22)/sx, 
y2  =  -x22  -  x21*(xl2*ul  -  x21*x22)/d, 
y3  =  (x!2*ul  -  x21-x22)/d, 
y4  =  1. 

Note  that  y2  =  -x22  -  x21*y3.  Let  sy  =  -sign(y2)*v/l+y22+y32 
We  seek  the  second  reflection  H2  =  I  -  wT/g  so  that 


1 1 


There  is  no  need  to  change  the  top  row.  Proceed  as  before  to  obtain 
vl  =  0 

if  -y2/sy  <  0.5,  then  v2  =  sy  +  y2;  else  v2  =  (i+y32)/(sy-y2) 
v3  =  y3 
v4  =  1 
g  =  v2-sy. 

Remark.  Since  each  H  that  swaps  A^  and  can  be  represented  in 
the  form  of  (6)  (lemma  2),  it  is  worthwhile  to  see  what  the  reflection, 
or  product  of  reflections,  looks  like  in  this  form.  In  fact  we  make  use 
of  it  in  Section  4.5.  We  compute  the  corresponding  and  C2  in 
appendix  B  They  are  obtained  by  noting  that 


4.5.  Special  treatment  for  the  diagonal  blocks. 

From  (8)  the  new  diagonal  block  A  is  equal  to  W'A’V1  for  some  W 
Given  A,W  we  use  a  special  subroutine  (called  EQUDI,  see  Appendix  C) 

to  put  2x2  diagonal  block  A  =W*A*W"1  into  standard  form  and  effect 
the  associated  changes  in  the  corresponding  rows  and  columns  of  the 
Schur  form.  For  better  accuracy  we  derive  here  the  analytic 
formulae  for  the  transformation,  based  on  W  and  A  Recall  that 
all  =  a22. 

d  =  det(W), 

2  =  a21*wl2*w22  -  al2*wll*w21; 
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all  ♦  2/d 

a21-w222  -  al2-w212 
d 

Apply  (9)  in  Section  4.1  to  A=WAW-1  above  to  find 

u  =  al2*(wll-w21)*(wll+w21)  +  a21,(w22-wl2Mw22+wl2) 

6  =  J4‘tan_1(2z/u) 
s  =  sin(d),  c  =  cos<6),  and 

all  vl 

PAP  = 

v2  all 

where 

-c  s 
s  c 

vl  =  [  a21*w22*(w22-wl2*c/s)  -  al2,w21,(w21-wll-c/s)l/d, 
v2  =  l-a21*w22*(w22+wl2*s/c)  ♦  al2*w21*(w21+wll*s/c)]/d 

Since  eigenvalues  are  preserved  under  similar  transformation,  we 
must  have  vl*v2=al2*a21;  thus  we  may  recompute  vl  from  v2  or 
vice  versa,  depending  on  which  one  is  smaller  in  magnitude. 


5.  Numerical  Tests 

We  have  done  extensive  testing  on  matrices  with  various  mixtures  of 
block  size.  All  3  algorithms  perform  well  in  most  cases.  To  investigate 
more  closely  the  accuracy  of  Algorithms  0,1,  and  2  under  extreme 
conditions,  we  tested  them  on  three  sets  of  matrices:  one  with  huge 

B,  one  with  a  choice  of  B  so  that  |det(X)l  «  11X11^,  and  finally  one  with 
fairly  close  eigenvalues.  These  tests  perform  2x2  block  swaps 

How  to  measure  the  ‘correctness”  of  the  computed  output  is  not  so 

easy.  Let  A  be  the  output  matrix  PTAP  where  P  is  the  orthogonal 
matrix  that  accumulates  all  the  transformations  that  are  applied  to 
A.  We  believe  that  the  only  sensible  measures  for  the  accuracy  of  P 

and  A  are  1)  how  close  is  P*PT  to  the  identity  matrix?  and  2)  how 

close  is  PAP1  to  the  original  matrix  A 


WAW  1  = 
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Thus  our  measuring  parameters  are. 

1.  Ep  =  II  I  -  PPT||,  (12) 

2.  Ea  =  IIA  -  PAPtH/IIA!|.  (13) 

Ep  is  the  orthogonality  error  in  P;  Ea  is  the  norm  relative  error  in 

PtAP.  Out  of  curiosity  we  also  computed 

3.  c  =  Maxflcjjl,  A(iJ)*0)  (14) 

where  Cjj  =  |(A  -  PAPT)(i,j)/A(ij)|.  This  is  the  worst  relative  error 
among  the  elements  of  PAPT. 

The  third  parameters  make  sense  only  when  A(i,j)  *  0,  since  fill-in 
(zero  elements  become  non- zero)  is  unavoidable  in  recovering  A  from 
P  and  A.  We  should  point  out  that  c  is  too  exigent  a  measure  for  the 
accuracy  of  P  and  A.  It  is  unreasonable  to  demand  high  relative 

accuracy  for  tiny  elements  in  PAPT.  Nevertheless  we  found  e  helpful 
in  showing  subtle  differences  between  good  swapping  programs  The 
following  results  were  obtained  on  a  SUN  3/50.  P  and  A  are 
computed  solely  in  single  precision  arithmetic.  However,  the  error 
measures  are  computed  in  double  precision.  Only  three  digits  are 
displayed  for  the  error  measures  in  order  to  keep  the  display  clean. 
We  have  also  run  our  program  on  a  VAX/750  with  similar  results. 

The  Fortran  program  for  Algorithm  0  is  Stewart’s  EXCHNG,  the 
Fortran  program  for  Algorithm  1  and  2  are  written  according  to 
section  4.3  and  4.4  with  special  formula  for  the  diagonal  blocks 
described  in  section  4.5.  See  the  listing  in  Appendix  C. 

The  roundoff  unit  is  2"^3  £  1.192E-7  in  the  following  numerical 
results. 

Test  matrix  1  with  parameter  t  (large  B) 

2  -87  -20000t 

5  2  -20000t 

A(t)  = 

0  0  1 

0  0  37 


IOOOOt 

-IOOOOt 

-11 

1 


14 


T 

Algorithm 

e 

see  (14) 

Ea 

see  (13) 

EP 

see  (12) 

1 

0 

2.11e-3 

1.45e-6 

1.66e-6 

1 

2.56e-6 

1.62e-7 

1.96e-7 

2 

2.74e-6 

3.57e-7 

3.41e-7 

2‘3 

0 

125e-4 

4.89e-7 

1.05e-6 

1 

1.66e-6 

127e-7 

1.62e-7 

2 

8.70e-6 

4.15e-7 

4.79e-7 

2-6 

0 

5.18e-5 

4.69e-7 

6.97e-7 

1 

1.90e-6 

2.19e-7 

2.32e-7 

2 

8.34e-7 

1.93e-7 

2.14e-7 

P  and  A  =  PTAP  from  algorithm  2  when  A=A(1). 


X  = 


P  = 


-57387.61  6294.046 

-3106.521  -7298.501  ' 

-9.999731E-1  7.337808E-3  -1.655211E-5  7.307688E-6 

-7.337804E-3  -9.999732E-1  -1.610292E-5  -1.307002E-4 
-1.675308E-5  -1.423457E-5  9.999108E-1  -1.334777E-2 
6.125366E-6  -1.309519E-4  1334783E-2  9.999109E-1 

1.000000  -85.98243  20011.92  -10194.38 
4.733524  1.000000  19985.38  9807223 

0.000000  0.000000  2.000000  -11.01783 
0.000000  0.000000  39.48143  2.000000 


Test  matrix  II  with  parameter  t  (Idet(X)  <<  IIXII2) 


A(-r)  = 


-3 

-87 

3576t 

4888t 

5 

-3 

-88^ 

-1440t 

0 

0 

17 

-45 

0 

0 

37 

17 

15 


T 

Algorithm 

e 

see  (14) 

Ea 

see  (13) 

EP 

see  (12) 

1 

0 

3.32e-5 

1.74e-7 

4.42e-7 

1 

i.99e-3 

i.93e-6 

6.39e-6 

2 

1.52e-5 

2.18e-7 

2.10e-7 

2-3 

0 

1.37e-5 

2.65e-7 

4.94e-7 

1 

2.41e-4 

9.09e-7 

1.85e-6 

2 

1.50e-6 

i.96e-7 

5.12e-7 

2-6 

0 

7.54e-6 

6.16e-7 

5.69e-7 

1 

7.94e-7 

1.99e-7 

2.60e-7 

2 

1.68e-6 

2.12e-7 

7.49e-7 

P  and  A  =  PTAP  from  algorithm  2  when  A=A(1). 


X  = 


P  = 


91.79044  -1242202 

-9.375406  19.85028  ' 

-9.907388E-1  -1.318285E-1  3216357E-2  -4.818546E-3 
1.35628  IE- 1  -9.640296E-1  2282810E-1  1.181298E-2 

4287299E-3  1.861262E-1  8.120817E-1  -5.530479E-1 

-4.807638E-3  1.364735E-1  5.360752E-1  8.330517E-1 

17.00000  -1432.443  -5568.100  -2230.135 
1.162349  17.00000  91.36795  824.0479 

0.000000  0.000000  -3.000000  -236.8775 
0.000000  0.000000  1.836391  -3.000000 


Test  matrix  III  with  parameter  t  (close  eigenvalues) 


A(t)  = 


7.001  -87  39  4t  222t 

5  7.001  -122t  -36t 

0  0  7.01  -11.7567 


0 


0 


37 


7.01 


Algorithm 
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e  E^  Ep 

see  (14)  see  (13)  see  (12) 


1 

0 

1 

2 

5.85e-6 

3.50e-7 

4.69e-7 

8.73e-7 

2.91e-7 

1.19e-7 

8.34e-7 

2.67e-7 

2.40e-7 

2'3 

0 

6.12e-6 

7.73e-7 

1.30e-6 

1 

6.12e-7 

2.78e-7 

2.36e-7 

2 

6.18e-7 

2.68e-7 

8.03e-7 

2'6 

0 

4.61e-5 

6.49e-7 

7.89e-7 

1 

320e-6 

4.45e-7 

4.05e-7 

2 

2.83e-6 

3.45e-7 

5.01e-7 

P  and  A  =  P^AP  from  algorithm  2  when  A=  A(l). 


X  = 


P  = 


1258173  -2869.060 

1218.421  1770267 

-1.0000Q0E+0  -1.293420E-4  6.833661E-5  -4.892B22E-5 
1293048E-4  -9.999997E-1  1.150727E-4  5.055802E-4 

6.830178E-5  1.152873E-4  1.000000E*0  4.074871E-4 

-4.902146E-5  5.055270E-4  -4.076063E-4  9.999997E-1 


A  = 


7.010000 

-87.01575 

-39.38432 

-22.17753 

4.999070 

7.010000 

12.19859 

36.00098 

0.000000 

0.000000 

7.000999 

-11.75932 

0.000000 

0.000000 

36.99190 

7.000999  . 

6.  Conclusion 


The  test  results  in  section  5  reveal  that  all  three  aglorithms  are 
acceptable  since  Norm  error  measures  E^  are  tiny  Algorithm  1  and  2 
have  the  advantage  of  keeping  the  real  eigenvalues  on  the  diagonals 
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at  all  times.  The  finer  measure  c  indicates  that  Algorithm  0  and 
Algorithm  1  in  certain  cases  are  inferior  to  Algorithm  2  but  in  other 
tests  cases  the  roles  of  Algorithm  1  and  2  are  reversed. 

We  find  no  reason  to  reject  any  of  the  methods  and  can  give  no 
perference. 
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Appendix  A.  Solving  AjX  -  XA2  =  B 


When  A}  and  A?  are  in  standard  form  the  inverse  of  the  coefficient 
matrix  can  be  expressed  quite  succinctly  and  safely.  Let 


A  = 

i 


«  8 

i  i 

If  <x 

i  i 


'  i=l,2.  If  <  0. 


Let  8=ai"«2.  then  the  equations  for  solving  X  may  be  written  as 


(1) 


c  \\ 
V2  c 


xll 

bll 

xl2 

bl2 

•x  =  b ;  x  = 

x21 

,  b  = 

b21 

x22 

b22 

where 


8  -* 

(2) 

C  = 

2 

K  6  J 

,  c2  = 

SH*2 

-2SB 


-28* 

2  2 

8 

2  2 


Multiply  (1)  as  indicated  in  order  to  make  the  coefficient  matrix  block 
diagonal. 


(3) 


°2-  Vt  ,  0 

•X  s 

C 

0  <r- 

V2  CJ 

Now  let 


G  =  (C2-  M  r1  = 
1  1 


t  26*, 

i 

288  t 

2 


/d 


where 


(4)  t  =  82+  M  -  6  1  .  d  =  t2-  (28p  )(28f  )  >  0, 
2  2  11  2  2 


and  premultiply  (3)  by  diag(G,G)  to  find 

<3  0 


x  = 


0  G 


ll  c 

-*  I  C 

J  1  12 
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where 


(5) 


G  0 
0  G 


ns 


'T\ 

-28v2 


♦*2  'T81  '2S<2P1 
1)6  "26p,  ?2  -T0j 

“av2  n8 

-**l  4*2  1)* 


n  =  t-2v2P2  =  &2-(Wi+P2*2>  >  0, 
4.  =  262  -  t  =  S2  ♦  (M1-M2)  • 


j 

i 


Inevitably  (5)  is  Cramer's  rule  and  d  =  det(A^eI  -  Ie  A2)  so  that  d=0 
if  and  only  a*  =  a2,  PiTfi  =  ^2 


Remark.  One  step  of  iteration  refinement  may  be  needed  if  the 
structure  matrix  is  ill-conditioned.  We  form  the  residual  matrix 
R=B-(AiX-XA2).  If  R  is  large  relative  to  B,  then  using  (5)  again  to 
solve  for  the  correction  matrix  E*  from  A^EX-EXA2  =  R  and  refine  X 
by  subtracting  Ex  from  X. 
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Appendix  B.  Representing  reflectors  in  form  (6)  of  section  3 
From  (6)  in  section  3  we  have 


Thus  to  represent  the  reflectors)  in  4.4  we  need  only  to  compute 
and  C2  using  the  formulae  above.  We  skip  the  details  of  algebric 
manipulations  and  give  the  results  below. 

case  lxl:  C*  =  (sx),  C2  =  (-sx),  where  sx  =  sign(xli)*yi+xll2 

case  1x2:  C\  =  (sx),  and 

-xll  1 

x12--L  -ill  .  <^C2)  =  sx. 

u3  u3 

where 

sx  =  -sign(xl2)*N/l+x  ll2+xl22  , 

if  -xl2/sx  i  0.5,  then  u3  =  xl2+sx;  else  u3  =  (l+xll2)/(sx-xl2). 


case  2x2: 

xll-SLl*  22  xl2__x21_ 

ul*v2  ul*v2 

x21  -  &  x22  -  — 

v2  v2 

-sx  0 

C  =  ,  det(C  )  =  det(C  )  =  sx*sy, 

2  [  yl  -sy  ]  1  2 
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where  sx,ul,yl,y3,y2,sy,v2  are  defined  by 
sx  =  sign(xll)'Vj/l+xll2+x212  , 

If  xll/sx  i  0.5,  then  ul  -  sx-xll;  else  ul=(l+x212)/(sx+xll), 
yl  =  -(xll*xl2+x21*x22)/sx 
y3  =  (xl2*ul-x21*x22)/(ul*sx), 
y2  =  -x22-x21*y3, 

sy  =  -sign(y2)-  Jl+yZ2+yZ2  , 

if  -y2/sy  <  0.5,  then  v2  =  sy*y2;  else  v2  =  (I*y32)/(sy-y2) 
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Appendix  C.  Listing  of  Fortran  Subroutines 


Subroutine  SWAPB  (Algorithm  1) 

Subroutine  SWAPB  (Algorithm  2) 

Subroutine  HOUSE  (used  in  Algorithm  2‘s  SWAPB) 
Subroutine  EQUDI 
Subroutine  TXMXT 


ooo  ooon  ooooooooooooooooooooooooo 


SUBROUTINE  SWAPB(T,P,N,  J1,N1,N2,NT,NP> 

REAL  T<NT,N),P(NP,N> 

INTEGER  NP,NT,N,  J1,N1,N2 

GIUEN  T  IN  SCHUR  PORN  SURPB  SWAPS  ADJACENT  DIRGOhRL  BLOCKS  T1 
AMD  T2  IN  MATRIX  T  BEGINNING  IN  AOW  J1  BV  ORTHOGONAL  SIMILARITY 
TRANSFORMATIONS  THAT  PRESERVES  THE  SCHUR  FORM  OF  T.  THE 
DIMENSION  OF  BLOCK  TJ  IS  N1  BV  Nt  AND  T2  IS  N2  BV  M2.  THE 
PARAMETERS  IN  THE  CALLING  SEQUENCE  ARE  < STARRED  PARAMETERS  ARE 
ALTERED  BY  THE  SUBROUTINE) 

*T  THE  MATRIX  UHOSE  BLOCKS  ARE  BEING  SUAPPED. 

*P  THE  AAAAY  INTO  WHICH  TVE  TRAMSFORMAT I ONS 
ARE  TO  BE  ACCUMULATED. 

N  TIC  OAOER  OF  THE  MATRIX  T. 

J1  THE  POSITION  OF  TIC  BLOCKS. 

N1  SIZE  OF  TIC  FIRST  BLOCK. 

M2  SIZE  OF  THE  SECOND  BLOCK. 

NT  TIC  FIRST  DIMENSION  OF  TIC  ARRAY  T. 

NP  THE  FIRST  OlfCNSION  OF  THE  ARRAY  P. 

HETHOO: 

ALGOR ITIfl  1  OF  “PROGRAMS  TO  SWAP  DIAGONAL  BLOCKS*  WITH 
SPECIAL  FORMULA  FOR  TIC  DIAGONAL  BLOCKS 

SUBPROGRAMS: 

TXMXT,  EQUDI 

INTERNAL  URRIABLES: 

REAL  D,R, S,Y,Z,U1,U2,U3, 01, 02,03, VI, V2, V3,Y4,W1,U2, W3,U4 

REAL  TI1,T22,T33 

REAL  X(2,2),X11,X12,X21,X22 

REAL  W(2,2>,M11,W12,W21,U22 

REAL  0<2,2>,011,012,021,022 

REAL  A1(2  2)  A2(2  2) 

EQUIOALENCE  «<  1,  1  >,X1 1 ),  <X<  1,2>,X12>,  <X<2,  1  >,X21  >,  <X(2,2>,X22  > 
EQUIOALENCE  <U< 1, 1 ),W1 1 ), <U( 1,2>,U12), CU(2, 1 ),U21 >, (W(2,2),U22 ) 
EQUIOALENCE  <0(1,  1  ),011  >,  <0<  },2),012>,  <0(2,  1  >,021  >,  <0(2, 2 >,022) 
INTEGER  IZ,  l,K,J1,J2,J3,J4 

SOLOE  X  FOR  II  -XI  IA1  01  II  XI  -  IA1  T121  BY  CALLING  TXMXT 

10  II  10  R21  10  11  10  R21 

CALL  TXNXT(T,N, J1,N1,N2,X, IZ,MT> 

IF  IZ-O,  A 1  AND  A2  ARE  TOO  CLOSE  TO  SWAP 

I F  < 12. EQ . 0  >  GOTO  50 
KH11H11+N2-2 
J2  =  J1+1 
J3  ■  JHN1 
04  ■  03+1 
IF(N1.EQ.2>  THEN 
Aid,  1>-T<J1,J1> 

A1(1,2KT<J1,J2> 

A1<2, 1 >=T(J2, J1 ) 

A1(2,2>-T<J2,J2> 


EM)  IF 

IF<M2.EQ.2>  THEM 
«2<1,1>»T<J3,J3> 

R2<1,2)=T<J3,J4> 

R2<2,1>T<J4.J3> 

R2<2,2)«T<J4,J4> 

ETC  IF 

GOTO  <10,20,30,40)  ,  K 
C 

C  H1«1,  N2=1  :  H  *  t  S  Oil  -X11  1  ),  S  *  1  O/SQRK 1+X1 1**2) 
C  f  0  S  ]  t  1  XI!  1 

C 

10  S»1.0/SQRT<1.0+X11*X11> 

C 

C  PERFORM  H*T*HT 
C 

Til  -  T<J1, Jt ) 

T22  «  T<J2,J2) 

DO  12  l*J1,N 
HI  «  T<J1, I ) 

U2  «  T<J2, 1 ) 

VI  -  S»<-X1 1*U1  ♦  U2 ) 

V2  *  S*<  HI  +  X11*H2) 

T<J1, I )  »  VI 
T<J2, 1  >  «  V2 
12  COMTIMUE 

DO  14  1  =  1,  J2 

HI  ■  T< I, J! > 

U2  =  T< I , J2> 

VI  -  S*<-X11*U1  ♦  U2 ) 

V2  =  S»<Ht  +  XI  1*W2) 

T<  I ,  J1 )  =  VI 
T< I , J2>  •  V2 
14  COMTIMUE 

C 

C  PERFORM  P=HT 
C 

DO  16  1=1, N 

HI  =  P<I,J1) 

U2  =  P< I , J2) 

vi  -  s*<-ximn  ♦  U2> 

V2  =  S*(  HI  +  xnno 
P<I,J1)  =  VI 
P<I,J2)  «  V2 
16  COMTIMUE 

C 

C  SURP  OIRGOMRL  ELEMENTS 
C 

T<J2, J2  >=T 1 1 
T(J1,J1)=T22 
GOTO  50 
C 

C  l  U1  0  011  -X11  1  0  ) 

C  M1«1,  N2»2  :  H  ■  1  U2  U3  Oil  -X12  0  1  1 

C  10  0  U1  1  I  1  XII  X 12 1 

C 

20  S  -  HX11*X11 


OOO  N  NOOU 


V  -  X  12*X  12 
U11«  SQRT<S> 

U22=  SQRT  < 1 . 0+V /S ) 

V21-  0 

U12-  XI  1*<X12AM  1  > 

U1  •  1.0/<SQRT<S+V>> 

Ul  ■  t.O/Ull 
U3  >  1.0/U22 
U2  -  -<X11*X12*U3>/S 
Til  *  T<J1,J1> 

J3  «  J2+1 

PERFORM  W*T*HT 

DO  22  l«U1,N 
Ul  «  T<J1, I ) 

U2  -  T<J2, I ) 

U3  *  T<J3, 1) 

VI  «  -X11*U1  ♦  U2 
V2  ■  -X12*U1  ♦  U3 
V3  «  Ul  ♦  XI  1*U2  ♦  X12*U3 
T<J1, I >  ■  U1*V1 
T<J2, I )  «  U2*Y1+U3*Y2 
T<J3, 1 )  *  U1*V3 
2  CONTINUE 

DO  24  l«1,J3 
Ul  «  T< I ,  J1 ) 

U2  «  T< I ,  J2> 

U3  =  T< I ,  J3) 

VI  -  -XI 1*U1  +  U2 
V2  «  -X  12*U  1  +  U3 
V3  *  Ul  ♦  XI 1*U2  ♦  X12*U3 
T<I,J1>  «  U1*V1 
T< I , J2)  =  U2*Y!+03*V2 
T< I , J3>  »  U1*V3 
4  CONTINUE 

PERFORM  P*HT 


DO  26  1  =  1,  N 

Ul  -  P<I,J1> 

U2  =  P<  I ,  J2 ) 

U3  =  P< I , J3) 

VI  ■  -XI  1*U1  ♦  U2 
V2  =  — X  12*14 1  +  U3 
V3  -  Ul  +  XI  1*U2  +  X12*U3 
PCI.J1)  =  U1*Y1 
P< I , J2)  «  U2*V1+U3*V2 
P< I ,  J3>  -  01*Y3 
26  CONTINUE 
C 

C  T<J3, J3>  ■  Til;  CfiLL  EQUOI  UITH  R2<2,2>,  U<2,2>  TO  GET  fl2' 
C 

T<J3, J3)  -  Til 

CflLL  EQUOI <T,P,N,J1,R2,U,U1 1*022, NT, NP ) 

GOTO  50 
C 


o  OO  <*»  W  OOO  OJOOOO 


I  U1  0  0  I  C  -X11  -X21  1  1 

M1*2,  N2*1  :  H  «  I  0  U1  011  1  0  Xlll 

(0  U2  U3  J  I  0  1  X21 1 

3  S  ■  1+X11*X11 
V  -  X2 1*X2 1 
Mil-  SQRT<S> 

U22»  SQRT( 1+V/S) 

M12-  0 

H21«  X11*<X21/U11) 

Ut  ■  1.0/<SQRT<S+V>> 

U1  -  1.0/U11 

U3  *  1.0/U22 

U2  -  -<X11*X2mi3>/S 

M2 1*  -U21 

D  *  U1 1 

Mil-  U22 

M22=  D 

T33=  T<J3,J3> 

PERFORM  W*T*HT 

00  32  l=J1,H 
U1  *  T(J1, I ) 

M2  -  T<J2, I ) 

M3  *  T<J3  I ) 

VI  *  -X 1 1*U 1  -  X21*M2  ♦  M3 
V2  *  HI  ♦  XI 1*U3 
V3  =  M2  +  X2t*W3 
T<J1, I )  -  U1*V1 
T<J2,I>  *  U1*Y2 
T(J3, I )  *  U2*V2+03*V3 
2  COHTIMUE 

DO  34  1*1. J3 
Ml  *  T< I , J1 ) 

M2  «  T< I ,  J2) 

M3  =  T< I . J3) 

VI  -  -X11*W1  -  X21*W2  ♦  M3 

V2  *  Ml  ♦  XI 1*M3 

V3  *  M2  +  X21*U3 

T< I , J1 >  -  U1*V1 

KI.J2)  *  U1*Y2 

T< I ,  J3>  *  U2*V2+U3*V3 

4  C0MT1MUE 

PERFORM  P*HT 


00  36  1  =  1,  M 

Ml  -  P<I,J1> 

M2  =  P<  1 ,  J2) 

U3  s  p(|  J3> 

VI  «  -X11*M1  -  X21*W2  ♦  M3 

V2  =  Ml  +  XI 1*M3 

V3  -  M2  ♦  X21*W3 

P<I,J1>  *  U1*V1 

P< I , J2 )  =  U1*V2 

P<  I ,  J3 )  -  U2*Y2HO*Y3 


36  CONTINUE 

C 

C  T<J1,J1>  ■  T33;  CALL  EQUDI  UITH  R1<2,2>,  U<2,2>  TO  GET  AT 
C 

T<J1,J1>  «  T33 

CALL  EQIDI<T,P,N,J2,A1,U,U11*W22,NT,1*>> 

GOTO  50 
C 

C  l  Ul  C  0  0  11  -III  -X21  1  0  1 

C  H 1*2,  N2*2  :  H  «  I  U2  U3  0  0  1  t  -X12  -X22  0  1  1 

C  f  0  0  Ul  0  1  t  1  0  XII  XI2  1 

C  [00U2U31I0  1X21X22) 

C 

40  CONTINUE 

0  «  X 1 1*X22-X 12*X2 1 
S  «  1+X1 1*X1 1 
0  -  X22*X22+0*0 
Z  «  X 12*X 12 
R  «  X21*X21 

V  »  S+Z 

U11  «  SORT(V> 

U22  -  SQRT<1.0+<D+fl)yV) 

U12  ■  0.0 

U21  «  <X 1 1*X2 1+X 12*X22 )/U  1 1 

V  -  S+R 

Ul  1  «  SQRT<V) 

U22  «  SORT <  1 . 0*<D*2  >/V  > 

V21  «  0.0 

U12  «  <X11*X12+X21*X22)AM1 
Ul  -  1.0  AM  1 
U3  ■  1.0AI22 
U2  «  -U12/<U11«U22> 

Ul  -  1.0/U11 
U3  *  1.0/U22 
U2  »  -U21/<U1ini22) 

U21  »  -M2 1 

V  «  Mil 
Mil  *  U22 
M22  *  V 

C 

C  PERFORM  H*mfT 

C 

DO  42  l=J1,M 

Ul  -  T<J1, I > 

U2  =  T(J2, I > 

U3  -  T<  J3,  I  ) 

U4  =  T<J4, 1  ) 

Vt  *  -Xtl*Ut  -  X21*M2  +  U3 
V2  -  -X12*U1  -  X22*U2  ♦  U4 

V3  ■  ui  +  xime  ♦  xi2nw 

V4  -  U2  ♦  X21*U3  ♦  X22*U4 
T<J1, I )  ■  U1*Y1 
T<J2,I>  ■  U2*V1  ♦  U3*Y2 
T< J3, I )  •  U1*V3 
T<J4, I )  «  U2*Y3  +  U3*Y4 
CONTINUE 
DO  44  1-1, J4 


42 


o  o  o 


Ml  -  T<I,J1> 

U2  *  T< I , J2> 

M3  «  T< I , J3> 
m  •  T< I , J4> 

VI  •  -X11*U1  -  X21*U2  ♦  W3 
V2  -  -X12mi  -  X22*U2  *  U4 
V3  *  mi  ♦  ximo  ♦  X12»W4 
V4  »  U2  ♦  X21*U3  ♦  X22*U4 
T< I ,  Jt  >  -  UI*V1 
T<I,J2>  »  U2*V1  ♦  U3*V2 
T<I,J3>  ■  U1*V3 
T<I,J4>  »  V2*V3  ♦  U3*Y4 
44  CONTINUE 

PEHFom  pnrr 

DO  46  l«1,N 

Ml  «  P<I,J1> 

M2  «  P< I , J2> 

M3  -  P< I ,  J3> 

M4  •  P<l , J4) 

VI  -  -X 1 1*M 1  -  X21*M2  +  M3 
V2  «  -X12*U1  -  X22*U2  *  M4 
V3  *  Ml  +  X11*U3  *  X12*W4 
V4  -  M2  ♦  X21*U3  +  X22*W4 
P<!,J1>  «  U1*V1 
P< I ,  J2>  *  U2*V 1  ♦  U3*V2 
P<I,J3>  «  U1*V3 
P< I , J4)  «  U2*Y3  ♦  U3*V4 
46  CONTINUE 

C 

C  CRLL  EQUOI  UITH  fll,  M  TO  GET  fir,  ft2,U  TO  GET  R2‘ 
C 

cfiu.  equoi <t,p,n, ji nt, w>> 

CflLL  EQUOKTjPjN,  J3,fl1,U,U1 1*W22,NT,NP> 

30  RETURN 

END 


SUBROUTINE  SUAPB<T,P,N,  J1,N1,M2,MT,NP> 

REAL  T<NT,N),P<NP,N> 

INTEGER  lf»,NT,N,J1,N1,N2 

C  GIUEN  T  IN  SCHUR  FORM  SWAPS  SWAPS  RDJRCENT  DlflGONRL  BLOCKS  T1 
C  AND  T2  IN  MATRIX  T  BEGINNING  IN  ROW  J1  BV  ORTHOGONAL  SIMILRRITV 
C  TRRN5F0H HT I ONS  THRT  PRESERVES  THE  SCHUR  FOFTI  OF  T.  THE 
C  DIMENSION  OF  BLOCK  T1  IS  N1  BV  N1  AND  T2  IS  N2  BV  N2.TFE 
C  PARAMETERS  IN  THE  CALLING  SEQUENCE  ARE  < STARRED  PARAMETERS  ARE 
C  ALTERED  BV  T)€  SUBROUTINE) 


c 

*T 

THE  MATRIX  WHOSE  BLOCKS  ARE  BEING  SWAPPED 

c 

*9 

TVE  AARRV  INTO  U4ICH  THE  TRANSFORMATIONS 

c 

ARE  TO  BE  ACCUMULATED. 

c 

N 

THE  ORDER  OF  THE  MATRIX  T. 

c 

J1 

THE  POSITION  OF  THE  BLOCKS. 

c 

N1 

SIZE  OF  THE  FIRST  BLOCK. 

c 

M2 

SIZE  OF  THE  SECOND  BLOCK. 

c 

NT 

THE  FIRST  DIMENSION  OF  THE  ARRAV  T. 

c 

A 

MP 

THE  FIRST  DIMENSION  OF  THE  ARRAV  P. 

c 

HETHOO: 

C  ALGORITHM  2  OF  "PROGRAMS  TO  SWAP  DIAGONAL  BLOCKS"  UITH 

C  SPECIAL  FORMULA  FOR  THE  DIAGONAL  BLOCKS 

C 

C  SUBPROGRAMS: 

C  TXMXT,  EQUDI 

C 

C  INTERNAL  UARIABLES: 

C 

REAL  X<2,2),U<4>,UU<4>,D,G,X11,X22,X12,X21,HRLF,V1,V2,V3 
REAL  U(2,2),U11,U12,U2t,U22 
REAL  U<2,2),U11,U12,U21,U22,TEIP,T11,T22,T33 
f£fS_  A  1(2  2)  R2<2  2) 

EQUIURLENCE  <X<1, 1  ).X1 1 ),  «<  1,2),X12>,  «<2, 1),X21),«<2,2).X22> 
EQUIUALENCE  <U<  1, 1  ),U11 ),  <U<  1,2>,U12>,  <U<2,  1  >,U21  >,  <U(2,2),U22> 
EQUIUALENCE  <V<1,  1  ),U11  ),  <U<  1,2>,U12>,  <U<2, 1  >,U21 ),  OK2,2),U22) 
INTEGER  K 
HALF  -  0.5 
C 

C  SOLUE  X  FOR  tl  -XI  IT1  T 12 1  (I  XI  =  IT  1  0  1  BV  CALLING  TXMXT 
C  10  II  10  T21  10  II  10  T2J 

C 

CALL  TXMXT<T,N,J1,N1,N2,X, 12, NT) 

C 

C  IF  IZ=0,  A 1  AND  A2  ARE  TOO  CLOSE  TO  SURP 
C 

IF< IZ.EQ.O)  GOTO  30 
K=Ht+N1+N2-2 
J2  -  J1+1 
J3  =  J1+NI 
J4  «  J3+1 
IF<N1.EQ.2>  THEN 
A1< 1, 1 >®T<J1, J1 ) 

A1<  1,2>»T<J1,  J2> 

AK2,1>T<J2,J1> 

A1<2,2)«T<J2,J2> 

END  IF 


00000  OOO  MOOO  ooo  -oooo 


IF<N2.EQ.2>  THEM 
R2<1,1>«T<J3,J3> 

R2<1,2)«T<J3,J4> 

R2<2,1)-T<J4,J3> 
ft2<2,2>»T<J4,J4> 

END  IF 

GOTO  < 10,20,30,40)  ,  K 

1.1  :  1*I-UU*/D,  H  t  XII  1  -  I  S  1,  S*S I GH<X 1 1 >*SQRT < 1+X 1 1**2 ) 
(II  (01 

3  S*  SIGM<SQRT< 1 .0+X1 1*X1 1 >,X1 1 > 

Til  *  T<J1,J1 ) 

T22  •  T<J2,J2> 

U<1)  «  S  -  XII 

IF«X11/$).GT.HflLF>  U<1>  «  1.0/<S+X11> 

U<2)  ■  1 
D  ■  U<1)*S 

CflLL  H0U$E<T,P,N, J1,U,2,D,NT,HP> 

SURP  DIRG0MRL  ELEMENTS 

T<J1,J1)  «  T22 
T<J2,J2>  *  Til 
GOTO  50 

1.2  :  I  1  XII  X12  1  H  *  l  0  0  S  1,  S»-SIGN<X12>*SqRT<  1+X1 1**2+X12**2> 

3  V  «  1  .CKX1 1*X1 1 

S  -  SIGN<SQRT<V+X12*X12>,-X12> 

U< 1 >  «  1 
U<2>  «  XII 
U<3>  -  X12  +  S 

IF<<-X12/S).GT.HflLF)  U<3)  »  V/<S  -  X12> 

D  »  U<3>+S 
VII  «  -X11 
V22  «  -XI 1/U<3) 

U21  -  1.0 

V12  »  -X12-1.0/IK3) 

Til  *  TCJ1, J1 ) 

CflLL  H0USE<T,P,H, J1,U,3,D,NT,NP) 

T<J3, J3)  ■  Til;  CflLL  EQU0I  UITH  ft2<2,2>,  U<2,2>  TO  GET  R2‘ 

T<J3,J3)  =  Til 

CflLL  EQUDl <T,P,N, J1,R2, V, S,NT,NP) 

GOTO  30 

2, 1  :  H  1-X111  -  I  SI,  S  -  S I GH<X 1 1 )*SQRT<  1+X 1 1*X1  HX2 1^X21 ) 

1-X211  *  l  01 
I  1  1  *  1  01 

30  T33  *  T(J3, J3> 

V  -  1.0+X21*X21 
S  «  SIGM<SQRT<V+XI1*X1 I >, X 1 1 > 

U<1)  *  S  -  XII 

IF<<X1 1/S>.GT.HflLF>  U<1>  -  V/<S+X11> 


OOOOO  £000000  ooo 


U<2>  -  -X21 
U<3>  «  1 
0  •  U<1  )*S 

U22  ■  X21AK1) 

U11  «  X21 

U12  -  -X11+1.0/U<1> 

U21  «  -1.0 

CfiLL  HQUSE(T,P,N, J1,U,3,D,NT,NP> 

T<J1,J1)  «  T33;  CfiLL  EQUDI  WITH  R1<2,2>,  W<2,2>  TO  GET  fir 
T<J1,J1>  »  T33 

CflU.  EQUOI<T,P,N, J2,R1,M,S,NT,MP) 

QOTO  50 

2,2  :  HI  1-X111  «  t  SI,  S  ■  SIGMCX1 1 >*$QRT< 1+X1 1*X1 1+X21*X21 > 

[-X211  -  I  01 

(  1  1  a  [  01 

t  0  1  a  (  01 

D  V  *  1 .  CHX2 1*X2 1 

SX  -  SIGM<SQRT<V+X1 1*X1 1  >,X1 1 ) 

U<1>  «  SX  -  XII 

IF<«11/SX>.GT.HflLF>  U<1>  »  V/<SX+X11> 

U<2)  -  -X21 
U<3>  *  1 
0  a  U<1)*SX 

CfiLL  H0USE<T,P,h,J1,U,3,0,HT,MP> 

TEMP  »  <T<J4,  J1  >*U<  1  HT<J4,  J2>*U<2HT(J4,  J3>*U(3>>/D 
T<J4,J1>  *  T<J4,J1>  -  TEMP*U< 1 ) 

T<J4,J2>  «  T<J4, J2)  -  TEMP*U<2> 

T<J4,J3)  a  T<J4,J3)  -  TEMP*U<3) 

VI  -  -<X1 1*X12+X21*X22)/SX 
V3  =  <X12m<1>-X21*X22)A) 

V2  a  -X22-X21*V3 

H2  IVIl  a  IV1 1  ,  HERE  V  a  H1*t-X121,  S  *  -SIGH<V2>*S0RT( 1+V2**2+V3**2 > 
IV21  -  f  SI  [-X22J 

IV31  »  [  oi  (01 

IV4 1  a  (  11  (11 

V  -  1 .0*V3*V3 
SY  a  SIGM<SQRT<V2*Y2+Y),-V2> 

ULK2)  *  SV+V2 

IF  <ABS<V2/SY > . OT , HflLF )  UU<2>aV/<SV-Y2> 

UU<3)  *  V3 
UU<4>  -  1 
G  «  UU(2>*SY 

CfiLL  HOUSE(T,P,H, J2,UU<2),3,G,MT,NP> 

TEMP  -  <UU<2>*T<J2, J1 HUUO >*T<J3,  J1  HUU<4)*T<J4<  J1  >>/G 
T<J2,J1)  a  T<J2,J1>  -  TETP*UU<2) 

T<J3, J1 )  a  T<J3,JO  -  TEMP*UUC3) 

T<J4,J1>  a  T<J4,J1)  -  TEMP*UU<4> 

U11  a  X22  -  1 .0/ULK2) 

W22  -  XII  -  <SV-X22>/<U<1>*UU<2>> 

U12  «  -X12+X21/<U<1>*UU<2>> 

U21  »  -X21+V3ALK2) 

VII  -  -SX 


V22  -  -SY 
V12  ■  Vt 
U2t  ■  0.0 
Y  -  SX*SY 
C 

C  CALL  EQUOI  WITH  A1<2,2>,  U<2,2>  TO  GET  fir,  A2,U  TO  GET  R2‘ 

C 

CALL  EQUDI<T,P,N, J1,A2,U,V,NT,NP> 

CALL  EQUOI <T,P,N,J3,A1,M,Y,NT,NP> 

50  RETURN 
END 

SUBROUTINE  HOUSE<T,P,N, J1,U,K,D,NT,NP> 

REAL  T<NT,N>,P<NP,N>,UOO,D 
INTEGER  K, J1,N,NT,HP 
C 

C  THIS  SUBROUTINE  PERFORMS  HOUSEHOLDER  TRRMSFOAMRT I  ON  ON  T  Att  ACCUMULATE 
C  THE  TRANSFORMATION  IN  P  : 

C  T  *  <I-UU*/D)T<I-UU*/D> 

C  P  *  P<l-W*/D>,  WHERE  U*  STANDS  FOR  THE  TRANSPOSE  OF  U. 

C  TIC  TRANSFORMATION  BEGINS  AT  T<J1,J1>.  THE  LENGTH  OF  U  IS  K. 

C 

C  INTERNAL  VARIABLES: 

C 

REAL  S, ZERO 
INTEGER  l,J 
ZERO=0.0 
C 

C  ROU  MODIFICATION 
C 

IF(O.EQ.ZERO)  GOTO  100 
00  30  J=J1,N 
S-0.0 

00  10  I  s  1 ,  K 

io  s«s+u<i>*t<ji+i-i,j> 

S-S/D 

DO  20  I-J1, J1+K-1 

20  T<  I ,  J  >T <  I ,  J  >-S*U<  I -J 1+ 1  > 

30  CONTINUE 

C 

C  COLUMN  MODIFICATION 
C 

DO  60  1*1, J1+K-1 
$-0.0 

DO  40  J*1,K 

40  S-S-HK  J  >*T  <  I ,  J 1+J- 1 ) 

S-S/D 

00  50  J*J1,  J1+K-1 

so  Ta^VTd^^snKj-Ji-M) 

60  CONTINUE 

C 

C  ACCUMULATION 
C 

DO  90  1-1, N 
S-0.0 

DO  70  J»1,K 

70  S-SHJ<J>*P<  I ,  J1+J-1 ) 


ooo  o o o  oooooooooonooo 


SUBROUTirt  EQU0KT,P,N,J1,fl,U,DETW,HT,lf>> 

PERL  T<NT,N>,P<NP,N),R<2,2>,U<2,2),0ETU 
INTEGER  J1,N,NT,NP 

THIS  SUBROUTINE  PERFORMS  R  UNITRRV  TRRNSFORMRT I  OH  <fl  REFLECTION) 
TO  NHKE  THE  DIRGONRL  ELEMENTS  IN  IWU*-!  EQURL  »C  PUT  THE  FINRL 
RESULT  IN  T. 

I-COS(Q)  SIN<Q>  1  l  Til  T12  1  t-COS<Q>  SIN<Q>  1  t  TIT  T12'  1 
I  SIN(Q)  COS(Q)  1  I  T21  T22  1  I  SIH<Q)  COS<Q)  1  =  t  T21*  TIT  J 
THE  TRRNSFORMRT I OH  IS  RCCUMULRTED  IN  < POSTMULT  I PLED  TO)  P.  THE 
INPUT  ft  MUST  HRUE  EQURL  DIROOHfiL  ELEMENT.  ICR E  OETW  IS  THE 
DETERM I NftNT  OF  U. 

F77  GENERIC  FUNCTIONS:  flTRN,  COS,  SIN 
INTERNRL  VARIABLES 
INTEGER  l , J, J2 

RERL  Q,U,S,C,Z,O1,02,TEMP,ZEBO,0r€,HflLF 
RERL  R11,R12,R21,R22,U11,U12,U21,U22 
fill  -  fl<1, 1) 

R12  «  fl<  1,2) 

R21  «  fl<2, 1) 

R22  -  fl<2,2> 

Ml  1  *  U(1,1) 

U12  *  U<  1,2) 

U21  ■  U<2, 1 ) 

1422  «  U<2,2) 

ONE  -  1.0 
HRLF  *  0.3 
ZERO  =  0.0 
J2-J1+1 

DETERMINE  THE  ANGLE  Q 

Z  =  R2mi12*W22  -  R12*U11*H21 

U  -  R12*<U1 1+1421  >*<U1 1-U21 )  ♦  R21*<U22-U12)*<U22+U12> 
IFCU+Z.EQ.Z)  THEN 
Q  *  RTRH< 1.0) 

ELSE 

Q  =  HRLF*flTHN<<Z+Z)/U) 

END  IF 

S  *  SIN<Q) 

C  «  COS<Q) 

NEW  DIRGONRL  BLOCK 

Z  -  C/S 

01  a  <  R21*U22*<W22-i!12*Z>-fl12*U21*<U21-il1  l*Z))/DETU 
Z  =  S/C 

02  «  <-R21'*W22*<U22+U12*Z>+A12,*H21*<U21+U1 1*Z))/DETU 
Z  *  R12*R21 

IF<RBS<01 ).GT.RBS<02))  TVCN 

02  «  z/vt 

ELSE 


01  -  Z/V2 


OOO  M  OOO  —  OOO 


END  IF 

ROU  NOOIFICRTION 

00  10  J*J2+1,N 

TEMP  -  -C*T<J1,JHS*T<J2,J> 
T<J2,J>-  S*T<J1,  JHC*T<J2,  J) 
T<J1,J)»  TENP 
0  CONTINUE 

COLDNN  NOOIFICRTION 

DO  20  1*1, J1-1 

TENP  -  T<I,J1>*<-C>4.T<I,J2>*S 
T<I,J2>*  TCI,J1)*S*TCItJ2>«C 
T<I,J1>®  TENP 
3  CONTINUE 

T<J1,J1>  «  fill 
T(J1,J2>  -  U1 
T<J2,J1>  •  U2 
T(J2.J2)  *  fill 

flccunuurriON 

DO  30  l»t,H 

TENP  =  P<I,J1>*<-C>+P<I,J2>*S 
P<I,J2>=  P<I,J1)*S+P<I,J2)*C 
P<I,J1>«  TE*f> 

30  CONTINUE 

40  RETURN 

END 


ooo  —  ooo  ooooonnoooooooooo 


SUBROUTINE  TXNXT<T,N,  IZ,NT> 

REAL  T<NT,N),X<2,2> 

INTEGER  N,NT,J1,N1,N2,IZ 

THIS  SUBROUTHC  SOLUES  FOR  HI  BV  N2  MATRIX  X  IM  T1*X  -  X*T2  =  T12. 
Tt  AND  T2  BEGIN  IN  ROUS  J1  AND  J2  RESPECT  I UELV,  T12  IS  THE  UPPER 
TRIANGULAR  PART  BETWEEN  T1  AND  T2.  THIS  PROGRAM  ASSUMES  THE 
DIAGONALS  OF  T1  <T2)  ARE  EQUAL.  THE  PARAMETERS  IN  THE  CALLING 
SEQUENCE  ARE  (STARRED  PARAMETERS  ARE  ALTERED  BV  TW  SUBROUTINE) 

*T  INPUT  MATRIX 

N  THE  ORDER  OF  THE  MATRIX  T 

C  J1  THE  POSITION  OF  THE  BLOCKS. 

N1  SIZE  OF  THE  FIRST  BLOCK. 

N2  SIZE  OF  THE  SECOND  BLOCK. 

•X  OUTPUT  MATRIX 

*IZ  OUTPUT  INDICATOR:  1-X  SOLUED  SUCCESSFULLY 

O-OUERFLOU  MRV  OCCUR. 

NT  THE  FIRST  0lf€NSI0N  OF  THE  ARRAY  T. 

INTERNAL  VARIABLES: 

REAL  D,DEL,BET1,BET2,GflM1,GAM2,T1,T2,TAU,E7R,PHI,DSC| 

REAL  P1,P2,P3,P4,P3,P6,P?,P8,P9 
INTEGER  I ,  J,K 
ZEAO-O.O 
IZ  «  1 

INITIALIZE 

DO  1  1-1,2 
DO  1  J*1,2 
X< I ,  J>=ZERO 
J2  -  J1+N1 

0EL»T<J1, J1  >-T<J2, J2) 

B11«T<J1,J2> 

IFCN1.EQ.2)  THEN 
J1P1  *  J1+1 
BET1  -  T(J1, J1P1 ) 

GAM1  »  TCJ1P1, J1 ) 

B21  =  T<J1P1,J2) 

END  IF 

IF<N2.EQ.2>  THEN 
J2P1  *  J2+ 1 
BET2  -  T<J2,U2P1) 

GAN2  =  T<J2P1,J2> 

BI2  -  T<J1,«J2P1 ) 

IF<N1.EQ.2)  B22*T <<J  IP  1 , J2P 1 ) 

END  IF 

K-M1+H1+N2-2 

GOTO  <10,20,30,40)  ,  K 

1  BV  1:  APH1*X  -  X*APH2  *  B1I 

10  IF(OEL.EQ.ZERO)  THEN 
IZ  «  0 

ELSE 


X<1, 1)  -  B1 1/DEL 


END  IF 
GOTO  50 
C 

C  I  BV  2:  APH1*IX1 1  X121  -  1X11  X121*lflPH2  BET21  ■  IB11  B12) 

C  [Gflft2  APH21 

C  _ 

20  0  ■  OEL*OEL  -  BET2^0FT12 

IFCD.EQ.2ER0)  TVCN 
12-0 

ELSE 

X< 1, 1 >  «  <OEL  *611  +  GflM2*B12>/D 
X<  1,2)  -  <BET2*B11  ♦  DEL  *612>/D 

END  IF 
GOTO  50 
C 

C  2  BV  1:  IflPHI  BET  1 1*1X1 1 1  -  IX111*HPH2  »  [Bill 
C  IQflNI  flPHI  1  tX21 1  1X211  tB21 1 

C 

30  D  =  0EL*0EL  -  BET1*GflM1 
IF<D.EQ.2£R0>  THEN 
12  *  0 

ELSE 

X<1,1>  »  <  DEL  <611  -  BET1*B21  )/D 
X<2, 1 )  =  <-GflH1*B11  +  DEL  *621 >/D 

EJCIF 
GOTO  50 
C 

C  2  BV  1:  IflPHI  BET  1 1*1X11  X121  -  1X11  X121*IflPH2  BET21  -  IB11  B121 
C  IGflni  flPHI 1  1X21  X221  1X21  X221  IGfl«2  APH21  IB21  B221 

C 

40  DSQ  «  DEL -DEL 
T 1  »  BET1*GflN1 
T2  -  BET2-0FTI2 
TflU  *  DSQ  +  <T2  -  T1 ) 

0  «  TAU*TflU  -  4*OSQ*T2 

IF(D.EQ.2ER0)  THEN 
12  *  0 

ELSE 

ETR  =  DSQ  -  <T 1+T2 ) 

PHI  *  DSQ  +  (T1-T2) 

T2  -  -<OEL+OEL) 

T 1  *  T2*BET1 
T2  *  T2*G»11 
PI  -  ETfl*OEL 
P2  »  PHI*GflH2 
P3  -  -TRU*BET  1 
P4  «  T1*GflH2 
P5  *  PHI*6ET2 
P6  -  T1*BET2 
P?  =  -TflU*GflM! 

P8  «  T2*Gflf12 
P9  -  T2*BET2 

X< 1, 1 >  «  <P1*ei1+P2*B12+P3*B21+P4*B22)/D 
X< 142>  -  <P5*611+P1*ei2+P6*B21+P3*B22>/D 
X<2, 1)  -  <P7*B11+P8*B12+P1*B21+P2*B22)/D 
X<2,2)  =  <P0*en+P7*ei2+PS*821+P1*B22)/D 
CC 


RBBB  RRR 


50 


COtFUTE  RESIDURL 

R1  »  Bit  -  <0EL*X<  1, 1  >-Gftt2*X<  t,2H6ETt*X<2, 1  >) 
R2  «  B12  -  <-BET2*X< 1, 1 HOEL*X< 1,2>+BET1^X<2#2>> 
R3  -  B21  -  <G«ni*X<1,  1H0EL*X<2,  1>-(3«12*X<2<2>> 
R4  •  B22  -  «3RM1*X<1,2>-6ET2*X<2, 1  H0EL*X<2,2>> 


PERFORM  OME  ITERRTIOM  IF  R*  IS  MOT  SMflLL  COMPRRED  TO  B* 

T1  «  RBS<B  11  HRBS<B12HRBS<B13HRBS<B  14  > 

T2  «  0.0625*<flBS<R1HRB$<R2HflBS<R3HflBS<R4>> 
IF<T1+T2.NE.T1)  THEM 

X<1,1>  *  X<  1, 1  >  ♦  <P1*R1+P2^2+P3^3+P4*R4)/D 
X<1,2>  *  X<1,2>  ♦  <P5*R  1+P 1*R2+P6*R3+P3*R4 )/D 
X<2,  1  >  -  X<2, 1  >  ♦  (PT^R  1+P(W2+P  1*R3*P2*fH  >/D 
X<2,2>  =  X<2,2>  +  <P9*R1+P7*R2+P5^3+P1*R4)/D 
END  IF 
EMDIF 
RETURN 
END 
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